Birth seasonality - Figures

1 Figure 1

official_birth_records = read_feather(path = str_c(IO$out_Rdata,"official_birth_records.feather"))

1.1 Births data

official_birth_records = official_birth_records %>% mutate(births_thousands = births/1000)

official_births_summary = official_birth_records %>% 
  group_by(country_area, country_area_col) %>% 
  summarize(min_births = min(births_thousands),max_births = max(births_thousands),mean_births = mean(births_thousands),
            max_date = max(date))
ref_date = as.Date("2020-01-01")

g_births = ggplot(official_birth_records, aes(x = date, y = births_thousands, col = country_area_col))
g_births = g_births +
  #### reference for variations
  geom_segment(data = official_births_summary, aes(x = ref_date, xend = ref_date, y = mean_births, yend = mean_births-5*mean_births/100),
               size = 1.5, alpha = 0.7)+
  # geom_text(data = official_births_summary, aes(x = ref_date, y = min_births-5*mean_births/100), label = "5% change from the mean",
  #           hjust = 0, vjust = 0, nudge_x = 100)+
  #### uncorrected births
  geom_line(aes(y = births_original_numbers/1000), col = "gray80")+
  #### corrected births
  geom_line()+   #geom_point(size = 0.5)+
  #### COUNTRIES
  geom_text(data = official_births_summary, aes(x = ref_date, y = max_births *1.05, label = country_area, fontface = 2), vjust = 0, hjust = 1)+
  #### settings
  scale_color_identity()+
  scale_x_date(date_minor_breaks = "1 year", date_labels = "%Y")+
  scale_y_continuous(expand = expansion(mult = c(0,0.2)))+ # sec.axis = sec_axis(~ (. - max(.))/max(.)*100, breaks = seq(-30,30,by = 10), name = "% change from the max") 
  ylab("Births (monthly in thousands)")+ xlab("Date")+
  facet_grid(country_area ~., scale = "free")+
  theme(strip.text.y = element_blank())

g_births

1.2 App data

### change
users = read_feather(path = paste0(IO$output_clue, "users.feather"))

n_users_per_country = users %>% group_by(country_area) %>% summarise(n_users = n()) %>% 
  mutate(country_area_col = dict$country_area$country_area_col[match(country_area, dict$country_area$country_area)],
         country_area = factor(country_area, levels = dict$country_area$country_area))


g_n_users = ggplot(n_users_per_country, aes(x = country_area, y = n_users/1000, fill = country_area_col))+
  coord_flip()+
  geom_bar(stat = "identity")+
  scale_fill_identity()+
  xlab("")+ylab("K# of app users")+
  theme(strip.text = element_blank(),
        axis.text.y = element_blank())+
  facet_grid(country_area ~ ., scale = "free")
g_n_users

1.3 Overall relative sex

pp = read_feather(path = str_c(IO$output_clue,"pop_indicators/relative_sexual_indicators_detrended.feather"))

pp_long = pp %>% filter(BC == "all") %>% 
  dplyr::select(date, country_area, n_users, n_sex, n_prot_sex, n_unprot_sex) %>% 
  rename(total_sex = n_sex, protected_sex = n_prot_sex, unprotected_sex = n_unprot_sex) %>% 
  pivot_longer(cols = c("total_sex","protected_sex","unprotected_sex"), names_to = "sex_type", values_to = "n_log") %>% 
  mutate(sex_type_str = sex_type %>% str_replace(.,"_"," ") %>% str_to_title())

pp_sum = pp_long %>% 
  group_by(date, sex_type, sex_type_str) %>% 
  summarize(n_users = sum(n_users),
            n_log = sum(n_log)) %>% 
  mutate(rel_counts = n_log / n_users) %>% 
  arrange(sex_type , date) %>% ungroup() %>% group_by(sex_type) %>% 
  mutate(
    t = row_number(),
    rel_trend = loess(rel_counts ~ t) %>%  predict(),
    sex_type_str = sex_type_str %>% factor(., levels = c("Total Sex","Protected Sex","Unprotected Sex")))


g_sex_freq = ggplot(pp_sum, aes(x = date, y = rel_counts, col = sex_type_str))+
  geom_line()+
  geom_line(aes(y = rel_trend))+
  ylab("Sex frequency (daily)")+xlab("Date")+
  expand_limits(y = 0)+  
  scale_color_manual(name = "",values = c("gray40","aquamarine3","indianred1"))+
  theme(legend.position = "top",
        legend.direction = "horizontal",
        legend.title = element_blank(),
        legend.background = element_rect(fill = "white", color = "transparent", size = 0.5)) #color = "gray90", 
g_sex_freq

pp_sum = pp_sum %>% 
  group_by(sex_type) %>% 
  mutate(
    mean_rel_counts = mean(rel_counts),
    detrended_counts = rel_counts/rel_trend * mean_rel_counts)

g_rel_sex = ggplot(pp_sum, aes(x = date, y = detrended_counts, col = sex_type_str))+
  geom_line()+
  ylab("Detrended sex frequency (daily)")+xlab("")+
  expand_limits(y = 0)+  
  scale_color_manual(name = "",values = c("gray40","aquamarine3","indianred1"))+
  theme(legend.position = "top",
        legend.direction = "horizontal",
        legend.title = element_blank(),
        legend.background = element_rect(fill = "white",color = "transparent", size = 0.5)) # gray90

g_rel_sex

1.4 Clue screen

# import png
clue_screen = image_read("../Figures Tables Media/Media/Clue_app_screens.png")

clue_screen_ratio = image_info(clue_screen)$height / image_info(clue_screen)$width 
padding = 20
g_clue_screen = ggplot()+ coord_fixed(ratio = clue_screen_ratio)+
  background_image(clue_screen)+theme(plot.margin = margin(t=padding, l=padding, r=padding, b=padding, unit = "pt"))

1.5 Assembled Figure 1

fig_1 = arrangeGrob(
  g_births,     
  g_clue_screen,   
  g_rel_sex,
  ncol = 2, nrow = 2, 
  widths = c(1.5, 1),
  heights = c(0.55,0.45),
  layout_matrix = cbind(c(1,1), c(3,4)))

fig_1 = as_ggplot(fig_1) + 
  draw_plot_label(label = letters[1:3], size = 15,
                  x = c(0, 0.6, 0.6), y = c(1, 1, 0.45))

fig_1

2 Figure 2

# load Figure 2 data

load(file = str_c(IO$p_outputs,"sex_models.Rdata"), verbose = TRUE)

## Loading objects:
##   sex_models
##   sex_models_df

plotlist = list()
n_items = c()
cas = unique(sex_models_df$country_area)
ok = foreach(ca = cas) %do% {
  cat(ca, "\n")
  
  # retrieving the model for this country and only plotting for BC all and sex_type all
  j = which((sex_models_df$country_area == ca) & (sex_models_df$BC == "all") & (sex_models_df$sex_type == "unprot_sex"))
  glm_sex_behavior = sex_models[[j]]$model
  
  # Visualisation of the coefficient
  g_coef = ggplot_sex_activity_glm_coefficient(model = glm_sex_behavior, show_intercept = FALSE, show_weekly_patterns = FALSE, ylim = c(-0.4, 1))
  g_coef = g_coef + 
    ggtitle(ca)+
    theme(axis.text.x = element_text(size = 0.8))
  print(g_coef)
  
  plotlist[[ca]] = g_coef
  n_items = c(n_items, nrow(g_coef$data))
}

## Brazil - Central-West

## Brazil - Northeast

## France

## United Kingdom

## United States - California

## United States - Northeast

fig_2 = ggarrange(
  ggarrange(
    plotlist[[which(cas == "Brazil - Central-West")]], 
    plotlist[[which(cas == "United States - Northeast")]], 
    labels = c("a","b"),
    ncol = 2, nrow = 1,
    widths = n_items[which(cas %in% c("Brazil - Central-West","United States - Northeast"))] / 
      sum(n_items[which(cas %in% c("Brazil - Central-West","United States - Northeast"))])
  ),
  ggarrange(
    plotlist[[which(cas == "France")]], 
    plotlist[[which(cas == "United Kingdom")]], 
    labels = c("c","d"),
    ncol = 2, nrow = 1,
    widths = n_items[which(cas %in% c("France","United Kingdom"))] / 
      sum(n_items[which(cas %in% c("France","United Kingdom"))])
  ),
  nrow = 2, ncol = 1
)

fig_2

3 Figure 3

G_table = read_feather(path = str_c(IO$out_Rdata,"Gestation_par_table.feather"))

# t = seq(0,4*pi, by = pi/20)
# f = sin(t)
# plot(t, f, type = "l")

t = seq(30*7, 45*7)

d = foreach(ca = G_table$country_area, .combine = bind_rows) %do% {
  this_ca = G_table %>% filter(country_area == ca)
  res = this_ca[rep(1,length(t)),]
  res = res %>% mutate(t = t, d = dnorm(t, mean = this_ca$G*7, sd = Gsd0), G = this_ca$G*7)
  res
}

br = seq(min(t), max(t), by = 14)

d = d %>% 
  mutate(country_area_wrapped = str_replace(country_area, " - ","\n"),
         country_area_col = dict$country_area$country_area_col[match(country_area, dict$country_area$country_area)])

g_d = ggplot(d, aes(x = t, y = d, fill = country_area_col))+
  geom_area(alpha = 0.6)+
  geom_vline(aes(xintercept = G, col = country_area_col))+
  scale_color_identity()+
  scale_fill_identity()+
  scale_x_continuous(breaks = br , labels = br/7)+
  scale_y_continuous(breaks = NULL)+
  xlab("G, gestation duration\n(weeks)")+
  ylab(expression(paste("Density, d",tau)))+  #expression(paste("Value is ", sigma,",", R^{2},'=0.6'))
  facet_grid(country_area_wrapped ~ ., scale = "free")+
  theme(strip.text.y = element_text(angle = 0, hjust = 0))
g_d

ggplot(d, aes(x = t, y = country_area, alpha = d))+
  geom_tile()+
  scale_y_discrete(limits = rev(levels(G_table$country_area)))+
  xlab("G, gestation duration\n(weeks)")+
  ylab("Density")+
  guides(alpha = FALSE)+
  scale_alpha_continuous(range = c(0,1))+
  scale_x_continuous(breaks = br , labels = br/7)

model = image_read("../Figures Tables Media/Media/Figure 3-01.png")
model_ratio =  image_info(model)$height / image_info(model)$width 
g_model = ggplot()+ coord_fixed(ratio = model_ratio)+
  background_image(model)



models = image_read("../Figures Tables Media/Media/Figure 3-02.png")
models_ratio = image_info(models)$height / image_info(models)$width
g_models = ggplot()+ coord_fixed(ratio = models_ratio)+
  background_image(models)

fig_3 = ggarrange(
  ggarrange(
    g_model, 
    g_d,
    labels = c("a","b"),
    ncol = 2, nrow = 1,
    widths = c(1,1)
  ),
  ggarrange(
    g_models,
    labels = c("c"),
    ncol = 1, nrow = 1,
    widths = 1
  ),
  nrow = 2, ncol = 1,
  heights = c(1,1)
)

fig_3

4 Figure 4

births = read_feather(path = str_c(IO$out_Rdata, "simulated_births.feather"))

births_STL = read_feather(path = str_c(IO$out_Rdata, "simulated_births_seasonal_trends.feather"))

opt_par_df = read_feather(path = str_c(IO$out_Rdata, "optimal_parameters_and_AIC.feather"))


ca_levels = dict$country_area$country_area
births = births %>% mutate(country_area = factor(country_area, levels = ca_levels))
births_STL = births_STL %>% mutate(country_area = factor(country_area, levels = ca_levels))
opt_par_df = opt_par_df %>% mutate(country_area = factor(country_area, levels = ca_levels))

4.1 time-series

births = births %>% 
  mutate(country_area_wrapped = str_replace(country_area," - ","\n") %>%  
           factor(., levels = dict$country_area$country_area %>% str_replace(.," - ","\n")))


g_ts = ggplot(births %>%  filter(sex_type == "unprot_sex", BC == "all"), aes(x = year_month, y = sim_births, col = model))
g_ts = g_ts +
  geom_line(aes(y = births), col = "black", size = 1)+
  geom_line(size = 0.5)+
  guides(col = FALSE)+
  xlab("Date")+ylab("")+
  facet_grid(country_area_wrapped ~ model, scale = "free" ,switch = "y")+
  theme(strip.placement = "outside",
        strip.text.y.left = element_text(angle = 0, hjust = 0.5),
        strip.background.y = element_rect(fill = "gray90", color = NA))+
  ggtitle("Measured vs. simulated births")

g_ts

4.2 seasonal trends

g_st = ggplot(births_STL, aes(x = month, y = sim_births_seasonal, col = model))
g_st = g_st +
  geom_hline(yintercept = 0, col = "gray80")+
  geom_line(aes(y = births_seasonal), col = "black", size = 1)+
  geom_line()+
  scale_x_continuous(breaks = seq(0,12,by = 3))+
  guides(col = FALSE)+
  xlab("Calendar months")+ylab("Seasonal trends")+
  facet_grid(country_area ~ model, scale = "free_y")+
  theme(strip.text.y = element_blank())+
  ggtitle("Seasonal trends")
g_st

4.3 AIC

opt_par_df = opt_par_df %>% group_by(country_area, BC, sex_type) %>% 
  mutate(AIC_best_model = min(AIC),
         dAIC = AIC - AIC_best_model)

g_aic = ggplot(opt_par_df %>%  filter(BC == "all", sex_type == "unprot_sex"),
               aes(x = model, y = country_area, fill = dAIC))
g_aic = g_aic +
  geom_tile()+
  scale_fill_gradient(name = "", low = hsv(0.58,0.4,1), high = hsv(0.02,0.9,0.9))+
  xlab("")+ylab("")+
  #guides(fill = FALSE)+
  facet_grid(country_area ~ model, scale = "free")+
  ggtitle("Diff. with best AIC")+
  theme(axis.text = element_blank(),
        strip.text.y = element_blank())

g_aic

g_best_aic = ggplot(opt_par_df %>%  filter(BC == "all", sex_type == "unprot_sex", model == "A"),
               aes(x = country_area, y = AIC_best_model))
g_best_aic = g_best_aic +
  geom_bar(stat = "identity")+ coord_flip()+
  xlab("")+
  facet_grid(country_area ~ model, scale = "free")+
  theme(axis.text = element_blank(), strip.text.x = element_blank())
g_best_aic

4.4 Phase and amplitude

opt_par_df = opt_par_df %>%  
  mutate(country_area_col = country_area_col %>%  factor(., levels = dict$country_area$country_area_col))

g_F = ggplot(opt_par_df %>% filter(model != "A", sex_type == "unprot_sex", BC == "all"), 
             aes(x = Tp, y = alpha, col = country_area_col))
g_F = g_F +
  geom_segment(aes(xend = Tp, yend = 0, linetype = model))+
  geom_point()+
  scale_color_identity(guide = "legend",
                       labels = dict$country_area$country_area, 
                       name = "countries/areas" )+
  scale_x_continuous(limits = c(0,1),breaks = seq(0,3/4,by = 1/4), labels = c("Jan","Apr","Jul","Oct"))+
  ylab("Relative amplitude")+ xlab("Fertility peak time")+
  coord_polar(theta = "x", start = 0)
g_F

opt_par_df = opt_par_df %>%
  mutate(
    country = str_split_fixed(country_area, " - ",2)[,1],
    Tp_seas = (Tp-10/365 - 0.5*(country == "Brazil"))%%1)


g_F_seas = ggplot(opt_par_df %>% filter(model != "A", sex_type == "unprot_sex", BC == "all"), 
             aes(x = Tp_seas, y = alpha, col = country_area_col))
g_F_seas = g_F_seas +
  geom_segment(aes(xend = Tp_seas, yend = 0, linetype = model))+
  geom_point()+
  scale_color_identity(guide = "legend",
                       labels = dict$country_area$country_area, 
                       name = "countries/areas" )+
  scale_x_continuous(limits = c(0,1),breaks = seq(0,3/4,by = 1/4), labels = c("Winter\nSolstice","Spring\nEquinox","Summer\nSolstice","Autumn\nEquinox"))+
  ylab("Relative amplitude")+ xlab("Fertility peak time")+
  coord_polar(theta = "x", start = 0)
g_F_seas

4.5 Assembling Figure 4 together

fig_4 = ggarrange(
  ggarrange(
    g_ts,
    g_st,
    g_aic,
    labels = letters[1:3],
    ncol = 3, nrow = 1,
    widths = c(3,1,0.7)
  ),
  ggarrange(
    g_F,
    g_F_seas,
    ggplot(),
    labels = c("e","f",""),
    ncol = 3, nrow = 1,
    widths = c(1,1,1), heights = 1
  ),
  nrow = 2, ncol = 1,
  heights = c(2,1)
)

fig_4

5 Saving to pdf

Figures_path = "../Figures Tables Media/Figures/"
scale = 2

ggsave(plot = fig_1, filename = str_c(Figures_path, "F1.pdf"), 
       width = 17.5, height = 12, units = "cm", scale = scale)

ggsave(plot = fig_2, filename = str_c(Figures_path, "F2.pdf"), 
       width = 17.5, height = 10, units = "cm", scale = scale)


ggsave(plot = fig_3, filename = str_c(Figures_path, "F3.pdf"), 
       width = 8, height = 7, units = "cm", scale = scale)

ggsave(plot = fig_4, filename = str_c(Figures_path, "F4.pdf"), 
       width = 17.5, height = 13, units = "cm", scale = scale)